1. Field of the Invention
The present invention relates to an image stitcher, an image stitching method, a computer readable recording medium having an image stitch processing program recorded thereon, for automatically joining from a plurality of microscopic images, for example, the images and acquiring a composite image having a high resolution and having a wide field of view.
2. Description of the Prior Art
[1] Description of Conventional Optical Flow Calculation Method
A technique for calculating an optical flow from two images, and registering the two images on the basis of the obtained optical flow has been known. Description is made of a conventional method of calculating an optical flow.
(1) Lucas-Kanade Method
A large number of methods of calculating an apparent optical flow of a moving object in a moving image have been conventionally proposed. The Lucas-Kanade method which is a local gradient method out of the methods is one of the best methods. The reason for this is that the speed of processing is high, implementing is easy, and the result has confidence.
As for the details of the Lucas-Kanade method, see an article: B. Lucas and T. Kanade, “An Iterative Image Registration Technique with an Application to Stereo Vision,” In Seventh International Joint Conference on Artificial Intelligence (IJCAI-81), pp.674-979, 1981.
The outline of the Lucas-Kanade method will be described.
When a gray scale pattern I=(x, y, t) of image coordinates p=(x, y) at time t is moved to coordinates (x+δx, y+δy) with its gradation distribution kept constant after a very short time period (δt), the following optical flow constraint equation (1) holds:
                                                                        ∂                I                                            ∂                x                                      ⁢                                                  ⁢                          δx                              δ                ⁢                                                                  ⁢                t                                              +                                                    ∂                I                                            ∂                y                                      ⁢                                                  ⁢                          δy                              δ                ⁢                                                                  ⁢                t                                              +                                    ∂              I                                      ∂              t                                      =        0                            (        1        )            
In order to calculate an optical flow {v=(δx/δt, δy/δt)=(u, v)} in a two-dimensional image, another constraint equation is required because the number of unknown parameters is two. Lucas and Kanade have assumed that the optical flow is constant in a local region of an object.
Suppose the optical flow is constant within a local region ω on an image, for example. In this case, a square error E of a gray scale pattern which is desired to be minimized can be defined by the following equation (2) when the following substitutions, are made:I0(P)=I(x,y,t),I1(p+v)=I(x+u,y+v,t+δt)
                    E        =                              ∑            ω                    ⁢                                    (                              [                                                                            I                      1                                        ⁡                                          (                                              p                        +                        v                                            )                                                        -                                                            I                      0                                        ⁡                                          (                      p                      )                                                                      ]                            )                        2                                              (        2        )            
When v is very small, the terms of second and higher degrees in Taylor's dilation can be ignored. Accordingly, at a relationship expressed by the following equation (3) holds:I1(p+v)=I1(p)+g(p)v  (3)
where g (p) is a linear differential of I1 (p).
The error E is minimized when the derivative of E with respect to v is zero. Accordingly, a relationship expressed by the following equation (4) holds:
                                                        0              =                            ⁢                                                ∂                                      ∂                    v                                                  ⁢                E                                                                                        ≈                            ⁢                                                ∂                                      ∂                    v                                                  ⁢                                                      ∑                    ω                                    ⁢                                                            (                                              [                                                                                                            I                              1                                                        ⁡                                                          (                              p                              )                                                                                +                                                                                    g                              ⁡                                                              (                                p                                )                                                                                      ⁢                            v                                                    -                                                                                    I                              0                                                        ⁡                                                          (                              p                              )                                                                                                      ]                                            )                                        2                                                                                                                          =                            ⁢                                                ∑                  ω                                ⁢                                  2                  ⁢                                                            g                      ⁡                                              (                        p                        )                                                              ⁡                                          [                                                                                                    I                            1                                                    ⁡                                                      (                            p                            )                                                                          +                                                                              g                            ⁡                                                          (                              p                              )                                                                                ⁢                                                                                                          ⁢                          v                                                -                                                                              I                            0                                                    ⁡                                                      (                            p                            )                                                                                              ]                                                                                                                              (        4        )            
Therefore, the optical flow v is found by the following equation (5):
                    v        ≈                                            ∑              ω                        ⁢                                          g                ⁡                                  (                  p                  )                                            ⁡                              [                                                                            I                      0                                        ⁡                                          (                      p                      )                                                        -                                                            I                      1                                        ⁡                                          (                      p                      )                                                                      ]                                                                        ∑              ω                        ⁢                                          g                ⁡                                  (                  p                  )                                            2                                                          (        5        )            
Furthermore, the optical flow can be found with high precision by Newton-Raphson iteration, as expressed by the following equation (6):
                              v                      k            +            1                          =                              v            k                    +                                    ∑                                                g                  k                                ⁡                                  [                                                            I                      0                                        -                                          I                      I                      k                                                        ]                                                                    ∑                                                (                                      g                    k                                    )                                2                                                                        (        6        )            I1k=I1(p+vk),gk=g(p+vk),I0=I0(p)
(2) Hierarchical Estimation Method
The largest problem of the gradient methods, including the Lucas-Kanade method, is that they cannot be applied to a large motion because a good initial value is required. Therefore, a method of producing images respectively having resolutions which differ at several levels like a pyramid hierarchical structure to solve the problem has been conventionally proposed.
Images having resolutions which differ at several levels are first previously produced from each of two consecutive images. An approximate optical flow is then calculated between the images having the lowest resolution. A more precise optical flow is calculated between the images having resolutions which are higher by one level. The processing is successively repeated until the optical flow is calculated between the images having the highest resolution.
FIG. 4, FIG. 3, FIG. 2, and FIG. 1 respectively illustrate an original image, an image having a lower resolution than that of the original image shown in FIG. 4, an image having a lower resolution than that of the image having a low resolution shown in FIG. 3, and an image having a lower resolution than that of the image having a low resolution shown in FIG. 2. In FIGS. 1 to 4, S indicates one patch.
An optical flow is gradually found from the image shown in FIG. 1 (an image in a hierarchy 1), the image shown in FIG. 2 (an image in a hierarchy 2), the image shown in FIG. 3 (an image in a hierarchy 3), and the image shown in FIG. 4 (an image in a hierarchy 4) in this order. In FIGS. 1 to 4, an arrow indicates an optical flow vector found for each patch.
However, the problem is that in a real image, there are few regions containing sufficient texture, so that a reliable optical flow is not obtained.
[2] Description of Optical Flow Calculation Method Developed by Applicant
An optical flow calculation method developed by the applicant of the present invention presupposes hierarchical estimation for producing images having resolutions which differ at several levels like a pyramid hierarchical structure to gradually calculate an optical flow. A method of calculating an optical flow conforms to a gradient method such as the Lucas-Kanade method. That is, it presupposes an optical flow estimation method using a hierarchically structured gradient method. The Lucas-Kanade method is used as the gradient method.
The optical flow estimation method developed by the applicant of the present invention is characterized in that an optical flow obtained in each of stages of the optical flow estimation method using the hierarchically structured Lucas-Kanade method is filled by dilation processing. This will be described in detail below.
One of advantages of the Lucas-Kanade method is that the result of tracking has confidence. Tomasi and Kanade have showed that the trackability of a region can be calculated from image derivatives as follows (C. Tomasi and Kanade, “Shape and Motion from Image Streams: a Factorization method-Part 3 Detection and Tracking of Point Features,” CMU-CS-91-132, Carnegie Mellon University, 1991).
From a 2×2 matrix of coefficients G in the following equation (7) having as elements the squares of differentials in the vertical and horizontal directions of an image ω in a region, its eigenvalues are calculated, thereby making it possible to determine the trackability of the region:
                    G        =                              ∑                          p              ∈              ω                                ⁢                                    g              ⁡                              (                p                )                                      ⁢                                          g                ⁡                                  (                  p                  )                                            T                                                          (        7        )            
When both the eigenvalues of the matrix G are large, the region changes in orthogonal directions, and can be uniquely positioned. Consequently, the confidence γ of the result of tracking can be obtained by the following equation (8) from the smaller eigenvalue λmin and a gray scale residual E between the regions after the tracking:
                    γ        =                              λ            min                    E                                    (        8        )            
The inventors and others of the present invention have developed a method of interpolating a region having low confidence using the result having high confidence in the same hierarchy in an optical flow. This uses the result in a hierarchy at the next coarser resolution level for only an initial value for tracking and does not utilize the result in a hierarchy at the current resolution level which is paid attention to. Alternatively, it is assumed that an optical flow in a region which is hardly textured has a value close to optical flows in its surrounding regions, to fill a flow field by morphology processing.
FIG. 5 shows how flow vector dilation processing is performed.
The left drawing represents a map of the confidence of a flow vector on a gray scale. It is assumed that the blacker the map is, the higher the confidence is.
Obtained flow is first subjected to threshold processing. A white portion is subjected to threshold processing because the result has low confidence.
The dilation processing of the result in the flow field is then performed as follows in imitation of hole filling processing using a morphology operation in a binary image. A flow vector u(i, j) in coordinates i,j of a region can be calculated, as expressed by the following equation (9), upon being weighted depending on confidence γ from flow vectors at its four vicinities:
                                                                        u                ⁡                                  (                                      i                    ,                    j                                    )                                            =                                                ∑                                      p                    ,                    q                                                  ⁢                                                                            γ                      ⁡                                              (                                                                              i                            +                            p                                                    ,                                                      j                            +                            q                                                                          )                                                              ×                                          u                      ⁡                                              (                                                                              i                            +                            p                                                    ,                                                      j                            +                            q                                                                          )                                                                                                  γ                    A                                                                                                                                                            (                                      p                    ,                    q                                    )                                =                                  (                                      0                    ,                    1                                    )                                            ,                              (                                  0                  ,                                      -                    1                                                  )                            ,                              (                                                      -                    1                                    ,                  0                                )                            ,                              (                                  1                  ,                  0                                )                                                                                                        γ                A                            =                                                ∑                                      p                    ,                    q                                                  ⁢                                  γ                  ⁡                                      (                                                                  i                        +                        p                                            ,                                              j                        +                        q                                                              )                                                                                                          (        9        )            
The processing is repeated until all regions having low confidence which have been subjected to threshold processing are filled. The filling processing is performed in each hierarchy. The flow vector u(i, j) in the coordinates i,j of the region may be calculated upon being weighted depending on confidence γ from flow vectors at its eight vicinities.
FIG. 6a illustrates an optical flow which has been subjected to threshold processing for an image in a hierarchy, and FIG. 6b illustrates an optical flow after filling. In FIG. 6a, arrows are optical flow vectors whose confidence is judged to be high by the threshold processing, and X marks indicate portions whose confidence is judged to be low.
[3] Description of Stitching Technique of Images
The applicant of the present invention has developed a technique (image stitching software) for calculating a geometric transform factor between images and stitching a plurality of images picked up on the basis of the calculated geometric transform factor.
In the image stitching software developed by the applicant of the present invention, when a plurality of images picked up in a two-dimensional arrangement are connected to one another to produce one composite image, the plurality of images picked up are arranged in a two-dimensional manner on a monitor screen such that the respective arrangement positions conform to an actual relative positional relationship. A pair of adjacent images from which an overlapped portion should be extracted is determined utilizing information related to the arrangement positions (the relative positional relationship between the images), to extract for the adjacent images the overlapped portion of the images. Utilizing the overlapped portion of the images extracted for the adjacent images, a correspondence between feature points is established for the adjacent pixels. The feature points between which a correspondence is established for the adjacent images are utilized, to calculate a geometric transform factor for the adjacent images and stitch all the images utilizing the geometric transform factor calculated for the adjacent images.
In the present conditions, when the plurality of images picked up are arranged in a two-dimensional manner on a monitor screen such that their respective arrangement positions conform to the actual relative positional relationship, a user sets the arrangement position for each of the images. Therefore, an operation for arranging the plurality of images picked up in a two-dimensional manner on the monitor screen is troublesome. Further, the user sets the arrangement positions of the images picked up one at a time, so that errors in the arrangement of the images are liable to occur.